Numerical Methods for a Model for Compressible Miscible Displacement in Porous Media

A nonlinear parabolic system is derived to describe compressible miscible displacement in a porous medium. The system is consistent with the usual model for incompressible miscible displacement. Two finite element procedures are introduced to approximate the concentration of one of the fluids and the pressure of the mixture. The concentration is treated by a Galerkin method in both procedures, while the pressure is treated by either a Galerkin method or by a parabolic mixed finite element method. Optimal order estimates in I? and essentially optimal order estimates in Lx are derived for the errors in the approximate solutions for both methods. Introduction. We shall consider the single-phase, miscible displacement of one compressible fluid by another in a porous medium under the assumptions that no volume change results from the mixing of the components and that a pressure-density relation exists for each component in a form that is independent of the mixing. These equations of state will imply that the fluids are in the liquid state. Our model will represent a direct generalization of the model [3], [4], [7] that has been treated extensively for incompressible miscible displacement. The reservoir ti will be taken to be of unit thickness and will be identified with a bounded domain in R2. We shall omit gravitational terms for simplicity of exposition; no significant mathematical questions arise when the lower order terms are included. Let c¡ denote the (volumetric) concentration of the /th component of the fluid mixture, i= 1,...,«. We assume that a density pi can be assigned to the /'th component that depends solely on the pressure p; moreover, we shall take this equation of state in the form (i.i) V = z<^ ri where z, is the "constant compressibility" factor [12, p. 10] for the /'th component. The assumption of miscibility of the components implies that the Darcy velocity of the fluid is given by (1.2) u = --Vp, M where k = k(x) is the permeability of the rock and p = p(c) = p(cx,...,cn) is the viscosity of the fluid. Assume that no volume change is induced by mixing the components and that a diffusion coefficient, which can combine the effects of Received April 14, 1982. 1980 Mathematics Subject Classification. Primary 65N30. 441 1983 American Mathematical Society 0025-5718/83 $1.00 + $.25 per page License or copyright restrictions may apply to redistribution; see http://www.ams.org/journal-terms-of-use 442 JIM DOUGLAS. JR. AND JEAN E. ROBERTS molecular diffusion and dispersion [13], exists that is component-independent; i.e., let (1.3) D = = <í>(jc) is the porosity of the rock. Then, conservation of mass of the /'th component in the mixture is expressed by the equation (1.4) *dr~+ v ' (c'p'") " v ■ (p where q is the external volumetric now rate, and c, is the concentration of the /'th component in the external flow; c, must be specified at points at which injection (i.e., q > 0) takes place, and c, is assumed to be equal to c¡ at production points. Carry out the differentiation indicated in (1.4), divide by p,, and use (1.1). The following equation results: (1.5) Vc,) -zlDvcl ■ Vp = c,q. If the components are of "slight compressibility" [12, pp. 10-11], then the term ZjCjU ■ Vp is effectively quadratic in the velocity, which is small in almost all of the domain, and can be neglected; we shall do so. The term -z^VCj ■ Vp = ztpk~]u ■ DvCj is small in comparison to the term u ■ Ve, that comes from the transport term, as both z, and D are small; thus, we shall neglect this term as well, so that we arrive at the equations (1.6) -£ + b(c)-£ + uvcV ■ (Dvc) = (c c)q. We shall assume that no flow occurs across the boundary: (a) u ■ v = 0 on dti, (I 13) v ' (b) (Dvccu) ■ v = 0 on8fl, where v is the outer normal to dti. In addition, the initial conditions (a) p(x,0)=po(x), xEti, ( ' ' (b) c(x,0) = c0(x), xEti, must be given. Note that b(c) = d(c) = 0 if the compressibilities z, and z2 vanish; thus, the model for the compressible problem converges to that of the incompressible model as the fluids tend toward incompressibility. The purpose of this paper is to formulate and analyze two numerical schemes for approximating the solution of the system (1.12)-(1.14). In both procedures the concentration equation (1.12b) is treated by a parabolic Galerkin procedure. In the first method the pressure equation (1.12a) is treated by a parabolic Galerkin procedure as well, and in the second it is treated by a parabolic mixed finite element technique. The analysis is given under a number of restrictions. The most important is that the solution is smooth; i.e., q is smoothly distributed, the coefficients are smooth, and the domain has at least the regularity required for a standard elliptic Neumann problem to have //2(ß)-regularity and more if the piecewise-polynomial spaces used in the finite element procedures have degree greater than one. We shall also consider License or copyright restrictions may apply to redistribution; see http://www.ams.org/journal-terms-of-use 444 JIM DOUGLAS, JR. AND JEAN E. ROBERTS only molecular diffusion, so that D = will be assumed bounded below positively as well as being smooth. Only the continuous-time versions of our methods will be studied here. The finite element methods will be formulated in Section 2. The scheme based on approximating the pressure by a Galerkin method will be analyzed in Section 3, and the one employing a parabolic mixed finite element method for the pressure will be treated in Section 4. The principal results of the paper are embodied in the L2-estimates for the error given in the inequalities (3.35) and (4.31) and in the L00-estimates (3.38), (3.39), and (4.32). The L2-estimates are of optimal order, and, except for a factor of log //"', so are the L°°-estimates. 2. Formulation of the Finite Element Procedures. Let h = (hc, hp), where hc and hp are positive. Let 911,, = 0 and with ch where q < 0. Equation (2.3) will be used to give the concentration in connection with both methods for approximating p and u. In addition, the initial approximate concentration ch(0) must be determined; several ways to specify ch(0) will be indicated after the requirements on c(0) — ch(0) become clear in the convergence analysis. The boundary condition was used in the derivation of the weak form of the equation. The Galerkin method for the pressure is based on the weak form of (1.12a) given by (2.4) (d(c)^j,v]j+(a(c)vp,Vv) = (q,v), v E Hx(ti), License or copyright restrictions may apply to redistribution; see http://www.ams.org/journal-terms-of-use COMPRESSIBLE MISCIBLE DISPLACEMENT IN POROUS MEDIA 445 for 0 < t < T. Let °Jih C lVx,0C(ti) be a piecewise-polynomial space of degree at least k associated with another quasi-regular polygonalization of ti; then, (2.5) mi\\v-vh\\Xtq, for v E Wk + x-q(ti). Then,/?,,: J -» %h will be determined as the solution of (2.6) (d(ch)^f,v)j+(a(ch)vph,Vv) = (q,v), vE%h,tEJ, starting from initial values ph(0) to be discussed later. The definition of the overall Galerkin procedure is completed by requiring that (2.7) u,, = -a(ch )vph. The parabolic mixed method that we shall employ as an alternative way to approximate the pressure is a simple generalization of a method introduced and analyzed by Johnson and Thomée [10]. Let V = [v E //(div; ti): v ■ v = 0 on dti) and W = L2(ti). Then, a saddle-point weak form of (1.12a) is given by the system (a) Id(c)-rr.w) + (V • u,w) = (q,w). wEW, (2.8) V di I (b) (a(c) ]u,v) (V ■ v, p) = 0, vEV. Now, let Vh X Wh be a Raviart-Thomas [14], [18] space of index at least â: associated with a quasi-regular triangulation or quadrilateralization (or a mixture of the two) of ti such that the elements have diameters bounded by hp. If ti is polygonal, then impose the boundary condition v ■ v = 0 on dti strongly on Vh. If dti is curvilinear and k = 1, use the boundary element described in Johnson-Thomee [10] and impose the boundary condition on Vh by interpolating the condition v • v = 0 at two Gauss points on each boundary edge. (This is a nonconforming feature; the analysis below does not cover this case explicitly, though it should extend without much difficulty to cover it.) Since the outer boundaries of petroleum reservoirs are not known exactly and since the regularity of the solutions of the differential equations is usually quite limited by the presence of sources and sinks, we shall assume ti to be polygonal and, for intuitive purposes, think of k as being either zero or one; thus, we shall, in particular, assume in the analysis that the boundary condition is represented exactly in Vh. The approximation properties for Vh X Wh are given by the inequalities (a) inf ||t; «A||0 = inf \\v üa||lW < M\\v\\k+Xhk+X, (2 9) (b) inf ||v • (v vh)\\0 < M{\\v\\k+X + ||v • v\\k+x}hk+x, »^Vh for v E V n Hk+X(ti)2 and, in addition for (2.9b), v • v E Hk+X(ti), and (2.10) inf \\w-wh\\0*ZM\\w\\k+xhk+x, wEHk+x(ti). License or copyright restrictions may apply to redistribution; see http://www.ams.org/journal-terms-of-use 446 JIM DOUGLAS, JR. AND JEAN E. ROBERTS The mixed method for the pressure equation becomes the finding of a map {uh, ph}: J -> VhX Wh such that (a) (rf(cA)^,w)+(v«A,w) = (?,w), wEWh, (b) (a(chyluh,v)-(v -v,ph) = 0, vEVh, for t EJ. Initial values must be specified for ph(0); consistent initial values uh(0) can then be computed from (2.1 lb). The Galerkin procedure defined by the combination of (2.3) and (2.6)-(2.7) represents an essentially traditional finite element approach to the nonlinear parabolic system (1.12), and it is a direct extension of a known method for treating the incompressible miscible problem [3], [4], [7], [8]. Thus, it is of some importance to establish that this, the simplest, finite element method for the compressible problem converges at an asymptotically optimal rate for smooth problems. It is well known that the physical transport dominates the diffusive effects in realistic examples of incompressible miscible displacement. In the liquid-liquid, compressible model studied in this paper, the transport will again dominate the entire process. Thus, it is more important to obtain good approximate velocities than it is to have extreme accuracy in the pressure. As in the incompressible problem [6], this motivates the use of the mixed method (2.11), now of parabolic type, in the calculation of the pressure and the velocity. The two procedures described above can easily be generalized to treat the «-component model. The pressure equation can be handled exactly as above, so long as the argument c in the coefficients a(c) and d(c) is interpreted as the vector {c,,. ..,c„_,}. Instead of a single concentration equation, there will be n — 1 such equations, and it is convenient to allow the approximation of c, to lie in GV)l„. Select 3w/3/ as the test function in (3.6) and observe the relation jt(a(ch)VTT.VTT) = 2^a(ch)v^,V^]j + ( |§(cj-|r Vit, Vît) = 2(a(ch)v,,V^) -(fc(ch)^v,,Vn) +(^(ch)^V,,V,). Thus, (3.7) ( ^í \ 1 Ida, ,3c \ + [{a(ch)-a(c)}vp,V^)+({d(ch)-d(c)}d-£,^ 1 /3a 2 a?' -M-SM'-S) License or copyright restrictions may apply to redistribution; see http://www.ams.org/journal-terms-of-use 448 JIM DOUGLAS. JR. AND JEAN E ROBERTS The only term in (3.7) that must be handled delicately is that involving a(ch) — a(c); it will be carried momentarily while the others are treated in a straightforward manner. We have assumed the coefficients d(c) and a(c) to be bounded above and below by positive constants independently of c. Thus, (3.8) d, dir dt K l_d_ 2 dt WVTT (a(ch)VTT, V>t) o lO.oc |V77||0 + dp di 0,x (llillo + HSllo) + + 9t) dt x 9w dt 9/ \Vtt\ O.oo + -[{a{ch)-a(c))Vp.V-^] Make the induction hypothesis that (3.9) HVir||0,x for some constant Kx. Apply (3.4), (3.5), and (3.9) to (3.8) to obtain the inequality Kt (3.10) 3/ + --(a(ch)VTT, Vit) K2\\\vnl + UWl + h2l+2 + h2k" \+e + 2[(a(ch)-a(c))vp,V dm "37 where K2 depends on Kx and ||3c/3?||0oo, ||3/>/3i||0,0 \\dp/dt\\k+x. Now, integrate (3.10) in time to see that (3.11) dJ'W^ 2dr+(a(ch)VTT,Vm)(t) Jo II "' o and